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Abstract: This study contains the derivation of an infinite space Green's 
function of the time-dependent radiative transfer equation in an anisotropi- 
cally scattering medium based on analytical approaches. The final solutions 
are analytical regarding the time variable and given by a superposition of 
real and complex exponential functions. The obtained expressions were 
successfully validated with Monte Carlo simulations. 
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1. Introduction 

The radiative transfer equation (RTE), an integro-partial differential equation, is widely used 
in many areas of physics to model the propagation of waves in random scattering media. Ex- 
amples are the light propagation in biological tissue [1], the neutron transport theory [2, 3], the 
electromagnetic waves propagation in plasmas or in the atmosphere. Although this equation is 
used since many decades complete analytical solutions for the steady-state and time domains 
are only available for the simplest cases [2-4]. Due to this fact the RTE is usually solved by 
numerical approaches such as the Monte Carlo method [5], the discrete ordinate method [6], 
the finite-element method [7] or is approximated by the diffusion equation [8]. In general, nu- 
merical approaches are in contrast to analytical methods computationally very demanding. 

Recently, elegant analytical approaches for solving the RTE in the steady-state domain were 
discussed in [9-11]. In principle, it is well-known that solutions obtained in the steady-state 
domain can be transformed into the frequency-domain resulting in a complex-valued absorp- 
tion coefficient. Then, by performing numerically an inverse Fourier transform of the solution 
in the frequency domain the solution of the time-dependent RTE is obtained. However, this 
procedure involves in many cases difficulties due to the fact that the Green's function of the 
steady-state RTE is already not a square integrable function [10]. Therefore the resulting time 
domain signal becomes a strongly oscillating curve which is inappropriate for an exact analy- 
sis or applications within inverse problems. Furthermore, the performance of the fast Fourier 
transform, which indeed would lead to decreasing calculation times, is circuitous if specified 
time values are required. In the publication of [12] the derivation of a cumulant solution for the 
time-dependent RTE was presented. 

In this article an infinite space Green's function of the time-dependent RTE is derived where 
the final solutions have an analytical dependence on the time variable and are well convergent 
due to the exponential decay of the time-dependent modes. Within the derivation no approxi- 
mations are made so that also the causality condition is satisfied; see the Results section. The 
solutions are successfully verified by comparisons with Monte Carlo simulations [1]. The de- 
rived Green's function cannot only be used for applications in infinitely extended scattering 
media but for all applications where the well-known approximated solution obtained from the 
diffusion equation is used [13], as was demonstrated for the case of isotropic scattering [14]. 

2. Theory 

The time-dependent radiative transfer equation for the radiance / = /(r, /i,f) caused by the 
internal source density Q = Q{r,t) in a spherical symmetric medium has the form 

where A is the scattering operator defined by 

A/(r,Al,f) = |/(s-s')/(r,M',Od's'. (2) 
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Here Gt = (J,, + CT.s is the total interaction coefficient, a„ is the absorption coefficient, is the 
scattering coefficient and = s f is the cosine of the angle between the direction of propagation 
§ and the unit vector of the position vector f . The function /(§ • §') describes a rotationally 
invariant phase function which can be expanded in Legendre polynomials Pi{x), where the 
expansion coefficients // are defined by 



f,=2nf_^f{li)P,{pi)Aii. (3) 



nl 
-1 

For finding the impulse response in space and time the source term becomes to 



Qi^^^)-^2m (4) 

In the following a plane wave decomposition of a rotationally invariant function in spherical 
polar coordinates (r, 0 , (j?) is derived. The radiance can be written as a three-dimensional Fourier 
integral 

/(r,cos0r,f) = j /(^,cos0k,f)exp(;k-r)d^^. (5) 

Note that the indices r and k of the angular variables 0^ and (p^ indicate the spatial and trans- 
formed space, respectively. By using the relation [15] 

,jkr.o.(6r-B^) ^ 4^ £ ^ y*^^^Q^^ (Pk)F/„,(0r, <Pr), (6) 

/=() m=-l 

where Yi,„{9,(p) are the spherical harmonics, and defining expansion coefficients 

Ii{k,t)=2n r I{k,cosek,t)Pi{cos9k)smekd9k (7) 
Jo 

yields the Legendre polynomial series 

^('"'^'0 = E^^'('-'0^/(m), (8) 
/=o "^"^ 

where the radial dependent kernel has the form of an inverse spherical Hankel transform 

Iii'-^t) = ri,ik,t)Mkr)k^dk. (9) 



27r2 Jo 

The belonging forward transform is found via orthogonality as 

I,{k,t)^47t{-j)' ri,{r,t)ji{krydr. (10) 

JQ 

The function ji{x) denotes the spherical Bessel function of the first kind. Inserting the series 
representation with the belonging radial dependent coefficients leads to a linear and time in- 
variant system (LTI system) 

+ ^JkI,-i{k,t) + aMk,t) + ^^jkl,+,{k,t) = 5{t)q,, (11) 

where Gi ~ Ga + {I — fi)Gs and qi = 8io. For the Henyey-Greenstein model [8] the expansion 
coefficients become fi =g' I G Nq. In P^y approximation (A^ odd) equations must be considered 
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only for / = 0, ...,A^ with the initial conditions /_i (A:,f) = In+i {k,t) = 0. The time-dependent 
impulse response of a LTI system can be obtained via the one-sided Laplace transform with 
zero initial conditions 



poo 

Ii{k,s)=^{Ii{k,t)}= / I,{k,t)e-"dt. 

Jo 



(12) 



Due to the above integral transforms the LTI system (11) is converted to a set of linear equations. 
After a few simple matrix conversions the system of linear equations becomes 



Q[Io{k,s)Mk,s),...,lN{k,s)]' = q, 



(13) 



where q = [1,0, ...,0]'^ and Q = (A + i/cl) • diag(VT, %/3, ...,V2/VTT). Here the matrix I 
denotes the identity matrix. The complex symmetric tridiagonal matrix A = A^ has the form 



(14) 
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The solution of the matrix equation Eq. (13) is formally obtained as [16] 



[k){k,s),h{k,s),...jN{k,s)]^ 



D 



I 



-1 



(15) 



Here D = diag ^ 



1 



The inverse Laplace transform can be performed using a 
well-known transform pair within the linear system theory for the so-called transition matrix 

exp(-Af) =^"H(iI + A)-'}, (16) 

Thus the time-dependent vector components become for f > 0 

[/o(^,f),/i(^,f),...,/iv(^,f)]^-cDexp(-cAf)q. (17) 

A simple but inefficient method for evaluation of the matrix exponential is given by the power 
series 

^2 



exp(— cAf ) 



fi^A-I-^^A- 



=0 



1! 



[cty 

2! 



A-- 



{cty 

3! 



A^±--- 



(18) 



which converge for an arbitrary square complex-valued matrix. Despite the fact that mathemat- 
ical software generally provide the matrix exponential it is maybe advantageous with respect to 
the computation time to avoid formally this function and performing the eigenvalue decompo- 
sition A = UAU^^ where U is the transformation matrix, which contains the eigenvectors |m„) 
(n = 0, 1, ...,A^) of A. The diagonal matrix A = diag(Ao,Ai, ...,AAf) (A„ e C) contains the eigen- 
values of the equation det(A — Al) = 0. Due to the above eigenvalue decomposition equation 
Eq. (17) simplifies to 



Mk,t)MKt),-MKt)Y =cDUexp(-cA/)U-'q, 



(19) 
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where the matrix exponential can now be evaluated as 



exp(— cAf) 



/ exp(-cAor) 0 ••• 0 

0 exp(-cAif) ••• 0 

\^ 0 0 ••• exp(-cAA^f) ) 



(20) 



Note that all expansion coefficients of the Legendre series Eq. (8) have a time decaying be- 
haviour which arises from the sign of the real part of the complex eigenvalues. After some 
algebraic conversions the time-dependent modes needed for the radiance Eq. (8) are obtained 
as 

c ^ 

^'(^'') = — T L («|Mo)(M«IOexp(-cA„f), (21) 
V 2/ + 1 „=o 

where the vector |mo) denotes the first column vector of the matrix U^'. Note that for obtaining 
an expression for the radiance caused by an isotropic emitting spherical shell source of the form 

dir-r') 

Qir,t) = \^5{t\ (22) 

where r' denotes the radius of the spherical shell, the source coefficients qi = 5/o in system 
(11) must be replaced by qi = {^j)' ii{kr')- For a source distribution with an arbitrary spatial 
dependence the resulting radiance can be obtained via convolution. 

Furthermore, the result for the fluence caused by an unidirectional source which illuminates in 
direction So is the same as for the obtained radiance caused by a isotropic source multiplied by 
the factor An. Then the above defined cosine of the angle between the direction of propagation 
and the vector of the position has now the meaning = So f . 

2.1. Numerical implementation 

For a fast and accurate evaluation of the obtained solution we replace the continuous Hankel 
transform over an infinite range, equation Eq. (10), by the corresponding finite version 

4,(f) =47r(-y)' rii{r,t)j,{^,kr)r^dr, (23) 
Jo 

where R equals the radius of a large sphere. Now the continuous wave number k in relation 
Eq. (10) becomes to discrete values i^^./ which are dependent on the transform order I and are 
given by the positive roots of the equation ji {t,kiR) = 0. Inverting the finite transform via the 
completeness relation of the spherical Bessel functions yields an alternative inverse transform 
for the integral representation Eq. (9) 

^'('•'0 = ^ 14/(0 (24) 

Due to the exponential decay of the time -dependent modes Iki{t) the solution, in general, con- 
verges already with a small number of discrete wave numbers. Additionally, we performed 
comparisons between the finite Hankel transform and the continuous version over an infinite 
range {R °o) which was evaluated by a Gaussian quadrature resulting in an excellent agree- 
ment. Note that especially for the evaluation of the fluence the required roots are easily given by 
the positive zeros of the sine function £,ito = kn/RV k G N. For every root the system of linear 
equations must be solved via eigenvalue decomposition of the matrix A. The whole procedure 
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can be easily implemented within a small MATLAB script, see e.g. on the internet [17]. In 
contrast to the steady-state theory here the eigenvalues and eigenvectors become additionally 
complex-valued. However, if an eigenvalue X„ becomes complex then X* is also an eigenvalue 
due to the symmetry to the real axis in the complex plane so that in all cases the functions 
//(r,f) are automatically real-valued. The fluence [5] which is given by the lowest order coeffi- 
cient /o(r, f) becomes the simple form 

<I>(r,f) - j Iir,^',t)dH' = ^ f W,o(f)sin(^,o'-). (25) 

Regarding the computation time the fluence in P25 approximation using 250 discrete wave num- 
bers can be evaluated for 500 time values in sa 1 s using a state of the art personal computer and 
a MATLAB script. Note that such an approximation order is much more than needed in most 
relevant cases. It is mentionable that the fluence in P3 approximation, which is a real alternative 
to the often used diffusion Green's function, can be evaluated completely analytically without 
a numerical eigenvalue decomposition. 

3. Results 

In this section the obtained solutions are vahdated against Monte Carlo simulations and partly 
the diffusion theory. Monte Carlo simulations are in the limit of an infinitely large number of 
simulated photons an exact solution of the RTE. The corresponding equations for the infinite 
space fluence and radiance within the diffusion approximation can be obtained from [8]. In the 
following comparisons a isotropic emitting point source distribution and the Henyey-Greenstein 
phase function are considered. The refractive index of the infinite medium is assumed as « = 
1.4. First, we compared the obtained expression for the fluence, equation Eq. (25), with the 
Monte Carlo method for an anisotropically scattering medium in Fig. 1. It can be seen that the 




derived solution and the Monte Carlo simulation converge to the same fluence. Due to the fact 
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that the two curves are practically the same, the relative differences between both methods are 
shown in the inset of Fig. 1. By increasing the number of simulated photons used in the Monte 
Carlo simulation these differences can be made arbitrarily small. For the next comparison 
the fluence for an higher absorbing medium and two different anisotropy factors is shown in 
Fig. 2. 




01 < ^ • • • • 1 

0 0.01 0.02 0.03 0.04 0.05 0.06 
t(ps) 

Fig. 2. Time-resolved fluence at r = 3.05 mm for two different anisotropy factors in a scat- 
tering medium with properties tTa = 0. 1 mm^ ' and CTj = 1 mm^ ^ . 

The derived fluence (solid curves) agrees again with the Monte Carlo method (symbols) 
for both anisotropy factors. Figure 3 shows the fluence obtained from the analytical solu- 
tion and the Monte Carlo method for a relatively small distance to the isotropic emitting point 
source. Additionally, the fluence within the diffusion approximation is included for two differ- 
ent diffusion coefficients D (D = l/(3c7i) and D = 1/{3g'^)) reported in literature [8], where 
<^s = ( 1 - ^) As and (Ji = jLi„ + ( 1 - g)ii,. 

The derived analytical Green's function of the fluence (solid curve) converges also for the 
relatively small source detector distances to the same result as obtained from the Monte Carlo 
simulation (symbols) whereas the diffusion approximation (dashed-curves) shows for both dif- 
fusion coefficients large differences compared to the radiative transfer calculations. 
Finally we give a comparison of the time-resolved radiance for three different directions of 
propagation. The results are shown in Fig. 4 for the angles 9 = 0.5°, 9 = 90.5° and 9 = 179.5° 
(indicated in the legend by 0 = 0, 7r/2 and n). 

Similar as for the fluence both methods result within the stochastic nature of the Monte Carlo 
simulations in the same radiance. 

4. Conclusion 

In this article time-dependent Green's functions of the RTE for the radiance and the fluence 
were derived and validated with the Monte Carlo method resulting in an excellent agreement. 
The obtained solution is completely analytical up to the P3 approximation and enables the 
consideration of an arbitrary phase function which must be expanded in the basis of Legendre 
polynomials. 

The relative differences between the independent solutions of the time-dependent RTE were 
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only caused by the statistics of the Monte Carlo simulations. The main computational cost of 
the derived Green's function is due to the numerical determination of several partial fraction 
coefficients which results in an eigenvalue decomposition of a complex symmetric tridiagonal 
matrix. 

The applications of the derived equations are manifold. They can be directly used for eval- 
uating measurements inside the investigated scattering media, like biological tissue or tissue 
phantoms [18, 19]. Further, the derived equations can be applied for obtaining solutions for 
boundary-value problems, e.g. the semi-infinite geometry, considering approximative bound- 
ary conditions [14]. In addition, the obtained solutions are important for validation of numer- 
ical methods for solving the RTE, e.g. Monte Carlo simulations. Finally, the derived Green's 
functions can be applied for the boundary element method to solve the transport equation for 
an arbitrary geometry [20]. 
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